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ABSTRACT 

We study the effects of diffusion on the non-linear corotation torque, or horseshoe drag, 
in the two-dimensional limit, focusing on low-mass planets for which the width of the 
horseshoe region is much smaller than the scale height of the disc. In the absence of 
diffusion, the non-linear corotation torque saturates, leaving only the Lindblad torque. 
Diffusion of heat and momentum can act to sustain the corotation torque. In the limit 
of very strong diffusion, the linear corotation torque is recovered. For the case of 
thermal diffusion, this limit corresponds to having a locally isothermal equation of 
state. We present some simple models that are able to capture the dependence of the 
torque on diffusive processes to within 20 % of the numerical simulations. 

Key words: planetary systems: formation - planets and satellites: formation. 



> 

o 
o 

> 

X 
C3 



1 INTRODUCTION 

Low-mass planets, embedded in circumstellar gas discs, are 



subject to orbital migration ( Goldreich & Tremaine 1979 



19801 through disc tides. Until they can acquire a massive 



gas envelope, planets comparable in mass to the Earth can 
therefore significantly move away from their place of birth. 
It was long thought that this Type I migration was always 



directed inward, and alarmingly fast (Ward 1997 Tanaka 
|et al.|2002 l, which suggested that the survival of long-period 
planets could be problematic within this theory. Recent 
work has indicated that the situation may be different if 
the disc's thermodynamics is taken into account properly 



( Paardekooper fe Mellema||2006a| |Kley & Cr ida 2008 |Kley 
et al.|2009| |, in which case outward migration is possible. It 
was subsequently realised that this migration behaviour is 



due to an entropy-related corotation torque (Paardekooper 
|fc Mellema||2008| |Baruteau fc Masset||2008al, and is non- 
linear in nature (iPaardekooper & Papaloizou 2008). 



In Paardekooper et al. (20101, hereafter Paper I, we 



studied the non-linear corotation torque, or horseshoe drag 
I Ward|1991 1, in the two-dimensional adiabatic limit. It was 
shown that this horseshoe drag is dominated by a term pro- 
portional to a radial entropy gradient in the unperturbed 
disc, and can lead to outward migration if the entropy de- 
creases outward. In the absence of any diffusion, however, 
the corotation torque is a transient phenomenon and prone 
to saturation dQuinn & Goodmanll 19861 |Ward||2007|). The 



necessary gradients (vortensity in the (locally) isothermal 
limit, entropy and vortensity for the non-isothermal case) 
within the horseshoe region need to be restored on a li- 
bration time scale in order for the corotation torque to 
be sustai ned (|Masset||2001| |Ogilvie fc Lubow|[2003[ ). It was 
shown in Paardekooper & Papaloizou ( 2008 I that for non- 



isothermal discs, both thermal diffusion and viscosity are 
needed to sustain the non-linear corotation torque. 

Apart from possibly keeping the corotation torque un- 
saturated, diffusion of heat and momentum can also, in the 
limit of high diffusion, act to reduce the non-linearity of the 



torque. It was already noted in Masset ( 2002 ) that there 
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exists a cut-off for the isothermal horseshoe drag at high 
viscosity. |Paardekooper fc Papaloizou| ( |2009a| ) showed that 
in this limit, the corotation torque returns to its linear value. 
For non-barotropic discs, the limit of high thermal diffusion 
effectively corresponds to the locally isothermal case. 

Depending on radial profiles of density, temperature, 
viscosity and thermal diffusion, the latter mainly determined 
by temperature and opacity, a wealth of Type I migration 
behaviour is possible therefore. While the sign of the unsat- 
urated torque is independent of planet mass, saturation does 
depend on the mass of the planet, which can lead to different 
migration directions for different planet masses. Moreover, 
since the level of saturation effectively determines whether 
the non linear corotation torque can dominate Type I mi- 
gration, it is important to study the effects of diffusion on 
low-mass planet migration. 

In this paper, we study Type I migration in non- 
barotropic discs in the presence of thermal and viscous dif- 
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fusion. We start in section [2] by reviewing our disc models 
and numerical methods. In section [3] we briefly review the 
results of Paper I on the unsaturated values of the torque. 
We then move on to the effects of viscous diffusion in section 
[4] by considering the simple isothermal case, first looking at 
saturation, and then moving to the limit of high diffusion 
in which the linear torque is recovered. We extend these re- 
sults to the non-barotropic case in section [5] We discuss our 
results in section |6] and conclude in section[7| 



2 NUMERICAL METHODS 

For the basic equations and details of the equilibrium disc 
models we refer the reader to Paper I. As in Paper I, we 
study two-dimensional power law discs in surface density 
£ (with index —a) and temperature T (with index —0). 
We study either isothermal discs, for which the pressure 
P = c 2 £, with c s the isothermal sound speed, or solve the 
energy equation and take a perfect gas equation of state, 
P — (7— l)e, where e is the internal energy density and 7 the 
adiabatic exponent. The scale height of the disc H = Cs/fl is 
taken to be small compared to the orbital radius. Typically 
h = H/r — 0.05 at the location of the planet. 

We have used the same numerical methods as in Paper I. 



set 



FARGO (Fast Advection in Rotating Gaseous Objects, Mas 
(2000a|bl) is based on the van Leer upwind algorithm on 



a staggered mesh. The time step calculation takes out the 
average angular velocity, which is dominated by Keplerian 
rotation, so that the time step is limited by the residual an- 
gular velocity rather than the much larger Keplerian veloc- 
ity. An energy equation solver was introduced in [Baruteau] 
& Massct (2008a). Diffusion in the energy equation is ob- 



tained by adding a term proportional to the Laplacian of 
the entropy: 



de 
dt 



+ V ■ ev = -(7 - l)eV ■ v + xeV log 



£"T 



(1) 



where v is the velocity and \ a thermal diffusion coefficient. 
This equation is strictly equivalent to an advection-diffusion 
equation for the gas entropy p/S 7 . Using this form of the 
energy equation is a common approach in studies of turbu- 
lent energy transport (see Shakura et al.||l978 l. We take x 
to be a constant throughout the disc. 

RODEO (ROe solver for Disc Embedded Objects, 



Paardekooper & Mellema (2006b I ) is a finite volume 
method, using an approximate Riemann solver ( Roe||1981 1 
to calculate fluxes between computational cells. Both codes 
were shown to yield similar values for the unsaturated horse- 
shoe drag in the adiabatic case (see Paper I), with RODEO 
usually producing slightly larger values (up to 10%). In 
RODEO, explicit heat diffusion is implemented by solving 



dT 



-V ■ (xVT). 



(2) 



For RODEO, we choose x to be a function of radius such that 
the initial temperature profile is a stationary solution. We 
have found that this only becomes important when the diffu- 
sion time scale over the horseshoe width becomes compara- 
ble to the dynamical time scale. Equation [2] is formally only 
correct for constant density, but therefore has the advantage 
of affecting the temperature only, while viscosity affects only 



density structures. This way, thermal and viscous diffusion 
can be treated separately. 

The two different approaches for obtaining diffusion 
were found to yield similar results, with differences less than 
10%. In both cases, the diffusion equations are solved ex- 
plicitly, which leads to a restriction of the time step. Note 
that we have not included viscous dissipation in the equa- 
tions: we assume that the background temperature profile is 
set by an equilibrium between viscous heating and radiative 
cooling ( Kley fc Crida|[2008 \ , and that viscous and thermal 
diffusion will act to restore the original profile. 

We take the viscosity v to be a power law in radius 
such that the initial surface density profile is a stationary 
solution. It is easy to verify that this requires v oc r" -1 ' 2 . 



3 UNSATURATED TORQUES 

Throughout this paper, we will work with a fixed gravita- 
tional smoothing length for the planet's potential, b/h = 0.4. 
We discuss possible effects of different smoothing lengths in 
section [5~7| In Paper I, we found the following unsaturated 
torques for this smoothing length, which we repeat here for 
convenience: the Lindblad torque 

7r L /r = -2.5 -1.7/3 + 0.1a, (3) 

the barotropic part of the horseshoe drag 

TlVbaro/ro = 1.1 ( I - a) , (4) 

and the entropy-related part of the horseshoe drag 

7r h s,cn t /T = 7.9^, (5) 

7 

where £ = /3— (7— l)a is the negative of the power law index 
of the entropy. The barotropic part of the linear corotation 
torque reads 



7r c ,lin,baro/ro = 0.7 [ - — (I 



(6) 



and the entropy-related part of the linear corotation torque 
is given by 



7rc,lin,cnt/r = 2.2 - — U 



1.4 



(7) 



As in Paper I, all torques are normalised to To = 
(g//i) 2 E p rpf2p, where subscripts X p indicates quantity X 
evaluated at the orbital radius of the planet r — r p . 



4 ISOTHERMAL DISCS 

4.1 A simple saturation model 



Masset (2001 1 discussed a model for corotation torque satu- 



ration. There are a few drawbacks to this approach, however, 
the major one being that the dependence of the corotation 
torque on background gradients has to be put in by hand. 
This actually leads to three possible saturation laws. An- 



other difficulty is that since Masset ( 2001 1 considers a con- 



stant viscosity v and a constant background surface density, 
the unperturbed model has a radial velocity profile (i.e. it is 
not in equilibrium). At large viscosities, the coorbital region 
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for low-mass planets can be deformed, resulting in a Type 
III flow structure. Below, we present a version of the |Masset| 
(2001 1 model that is tailored to the problem at hand. 



Assuming Keplerian rotation everywhere, the evolution 
of the surface density is given by 



dE _ 

dt ' 

|n-fi p 

2tt 



3 j^ 

r dr 






i d . 

(r 

rQ. dr v 

-e 1 n 



Vsn) ) + 



(8) 



where the second term on the right hand side accounts 
for changes due to the horseshoe turns (see |Masset||2001[ |. 
Primed quantities indicate they should be evaluated at 
r' = 2r p — r, i.e. on the other side of the horseshoe turn, 
and II denotes the rectangular function: 



n(t) = 



i |*| <§ 
o |t|>|. 



(9) 



The use of the rectangular function ensures that changes due 
to horseshoe turns only occur for \r — r p \/r p < x s , where x s 
is the non-dimensional half-width of the horseshoe region. 



Changing variables to y — E/fi and x — (r 
have 



>)Api 



dy(x,t) 



dt 



l«^( w (- M )-*M))n(£ 



l y(x,t) 



+ 



(10) 



Note that y is proportional to the inverse of the specific 
vorticity, or vortensity, under the approximation of Q being 
strictly the Keplerian angular velocity everywhere. Equation 
1 10| then shows that the vortensity changes with time due 
to diffusion (first term on the right-hand side) and due to 
advection in the horseshoe region (second term on the right- 
hand side). 

From now on, we will look for stationary solutions. 
Changing variables for a second time, with z = 2(l + x) 1 ' 2 — 
2, and / = (1 + z/2) 2a - 3 y(z), we find 



3^ P d 2 .f(z) 
r 2 dz 2 

30,, 



47T 
1 + 



/(" 



/(*; 



n 



(ii) 



where we have approximated the second term of equation 
10| for i< 1 (or, equiyalently, z <C 1). Here, z s = z(x s ). 



Masset 



Following 

and an odd part, /, given by 



/(*) = o (/W + /("*)) : 



and 



/(*)=^ (/(*)- /(-*))• 



(20011, we split / into an even part, /, 



(12) 



(13) 



We then obtain finally 

d 2 f(z) 



1 p a 'p | 
2tt^ d 



A /(*) + U - a */(*) n 



2;, 



= 0. (14) 



The even part of equation 14 gives / = / p , if we demand that 
/ is bounded for all z. Note that f — f p exactly corresponds 
to the unperturbed profile of y(x). Defining 



27T^ p 

the odd part of equation [14] gives 



d 2 f(z) 
dz 2 



''■U(z)+[--a)zf P ) =0, 



(15) 



(16) 



for ^ z ^ z s , subject to boundary conditions at z — and 

z = Zs- 

Since / is odd, we must have /(0) = 0. The most ob- 
vious outer boundary condition, to ensure continuity of /, 
is to take f(z s ) = 0. This way, vortensity perturbations are 
localised to the horseshoe region only. However, for x s <C h, 
we expect pressure effects to spread vortensity perturba- 
tions over a length scale H ( Casoli fc Masset|2009 \ . In other 
words, assuming Keplerian rotation everywhere means that 
vortensity perturbations will not be in pressure equilibrium, 
while for i s < fe we expect pressure equilibrium to be main- 
tained at all times. 

We can incorporate this in the simple model by demand- 
ing that: 



Q J 



rE or 



(17) 



However, using equation |17| in equation [8] makes the prob- 
lem very difficult to solve. We can still make progress, if we 
assume that Keplerian rotation still holds for x < x s , while 
we use equation fF7\ in equation [§] for x > x B . In order to 
match the two solutions at x = x s , we solve equation [16] 
with boundary conditions /(0) = 0, f(z s ) — f s , which gives 
us the inner solution: 



M*) = 




(18) 



where I denotes a modified Bessel function. We fix / s by 
demanding that the inner solution matches in a smooth way 
to the outer solution. From above equation, it is clear that 
saturation is governed by a single parameter p: 

p = 2x/fcr|/3. (19) 

The outer solution / ou t will decay on a scale H towards 



the equilibrium solution (Casoli & Masset 20091. There is 



no need to calculate it exactly, as long as we have a; s <C h. 
We know that the odd part of the derivative of / out will 
approach zero on a length scale H , so for x s ^C h we can set 
it to zero. We can therefore fix / s by demanding that the 
derivative of the inner solution is zero at x s . Note that this 



is the same boundary condition as used in Masset (20011 
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Figure 1. The function y(x) = T,(x)/il(x), for three different 
saturation parameters p. The inset shows a larger spatial domain, 
where the exponential decay of the perturbation is apparent. 



This leads to 
3 _ 
2 



f, 



s / p 



3P^4/3(P) 



(20) 



p 2/ 1/3 (p) + 3p/4/3(p)' 

The resulting solution for y(x) is shown in Fig. Ill for three 
different values of p, where we have chosen the outer profile 
for x > x s to approach the equilibrium profile exponen- 
tially on a scale H. The exact form does not matter since it 
does not affect the torque. Note that all profiles are indeed 
smooth at x — x s . 

The horseshoe drag in an isothermal disc is given by 

y^x 2 dx, 

Vp 
where y(x) is the odd part of y(x). 

The torque on the planet can be obtained by again in- 
tegrating over all material that has undergone a horseshoe 
turn. This is why the exact form of the outer solution is of 
no importance: since that material never executes a horse- 
shoe turn, it does not contribute to the torque. We then get 
a multiplier for the total horseshoe drag 

8/4/3 ip) 



-L hs,baro — OZ-fpT^L^ 



(21) 



F{p) = 



3 P /i/ 3 (p) + %V 2 h 



(22) 



H/s(p) 

Equation [22] is shown in Fig. [2] together with numerical re- 
sults obtained with FARGO. The model matches the numeri- 
cal results to within 20 % for p > 1. For smaller values of p, 
or larger values of z/ p , viscosity starts to affect the horseshoe 
turn itself, a regime we deal with in section |4.2| below. It 
can be shown that equation [22] is identical to equation 20 
in Masset (20021. Note, however, that in our case we auto- 



matically get the correct density profile dependence, which 



removes the ambiguity in the choice of F in Masset (2001 1 



We comment that F(p) may be approximated to within 5 % 
by a simpler form 

F(P) = T+W^F- (23) 



standard 
a p=l/2 
□ b/h=2 



Ph 




100.0 



Figure 2. The function F(p). The solid curve represents equation 
|22| Symbols indicate the results of numerical simulations using 
FARGO of a q = 1.26 • 10 -5 planet embedded in a h = 0.05 disc. 
The standard model has a = 1/2, /3 = 1 and b/h = 0.4. 



4.2 The limit of high v p 

In the limit p < 1, the assumption that the horseshoe 
turn over time is much smaller than the viscous diffusion 
time scale across the horseshoe region breaks down. In this 
regime, the horseshoe turn itself is modified by viscous ef- 
fects, and specific vorticity is no longer conserved during the 
turn. It was shown in |Paardekooper fc Papaloizou| ( |2009a[ ) 
that the non-linear horseshoe drag is then replaced by the 
linear corotation torque, which is in general much weaker. 
Because of the strong viscosity required, saturation does not 
play a role in this regime, making the torque fully linear. It 
is important to understand this regime, since in the non- 
barotropic case it corresponds to the transition between the 
optically thick, mostly adiabatic, inner disc, and the opti- 
cally thin, locally isothermal, outer disc, with viscosity re- 
placed by thermal diffusion. 

The regime in which the viscosity is strong enough to af- 
fect the horseshoe turns, but not strong enough to make the 
torque fully linear, is very difficult to model. All we can say 
is that in the limit p — > 0, there are basically no turns as far 
as the flow is concerned, since the viscosity is strong enough 
for the flow to forget where it came from, and the coro- 
tation torque will be fully linear. In the opposite limit, the 
horseshoe turns will prevent the linear torque from being set 
up, and the corotation torque will be given by the (possibly 
saturated) horseshoe drag. As the viscosity increases, the 
horseshoe drag decreases and the linear corotation torque 
increases, and it is not a-priori clear that these processes oc- 
cur at the same rate. In fact, one may expect the rates to be 
different, since the time scale to set up the horseshoe drag, 
a fraction of the libration time scale, is longer than the time 
scale associated with the linear corotation torque, which is 
a dynamical time scale. 

We therefore make the assumption that the corotation 



Non-isothermal Type I planetary migration -II 5 



torque is given by 

-L c,baro — tjl hs,baro i V / c,lin,baro. 



(24) 



with ^ G, K ^ 1. In the limit of low viscosity (but still 
high enough to keep the corotation torque unsaturated), we 
expect G,K — > 1, and the corotation torque is given by the 
horseshoe drag. In the limit of high viscosity, G, K — > 0, and 
the torque will be fully linear. As mentioned above, we allow 
the horseshoe drag to decrease at a different rate than the 
linear torque increases towards high viscosities by having K 
different from G. Below, we will derive a general form for G, 
from which we will extract K using a time scale argument. 
Consider the contribution of a single streamline, making 
a horseshoe turn from —x to x, to the horseshoe drag. It 
was shown in |Paardekooper fc Papaloizou| ( |2009a| ) that the 
time scale to set up the horseshoe drag is approximately 
Ttum = 3riib/20. When this is comparable to the diffusion 
time scale across x, viscosity is able to affect the horseshoe 
drag. We take this into account by taking the contribution 
of the streamline to the horseshoe drag to be 



dr h 



ro — O 



a ) ^r^fl^F ( t — 



d.i\ 



(25) 



with r = Tdiff /Ttum = A5irp 2 /4 (using r tu 



2/(5x s fip) and 



Tdiff = x B r p /u p ), and T{t) such that .F(O) = and J-{t) — > 1 
for i — Y oo. For rtum -C Tdiff, this gives us the ordinary 
horseshoe drag, while for 7t U rn 3> Tdiff, the contribution of 
the streamline goes to zero. 

To keep all integrals tractable, we choose a particularly 
simple form for T: 



T(t) = 



{t/roY t < to 

1 £ > To, 



(26) 



for some I, To > 0. We have obtained good agreement with 
numerical experiments with I = 3/4 and tq = 2. It is 
straightforward to show that in general 

G{t) = 4 / w s F{rw i )dw = ^t~ 4/3 I" t 1/3 T{t)dt, (27) 
Jo 3 J 

which for our particular form for T yields 

, i 



G(r) 



4+3! I T 



T < T 

This can be written in terms of p rather than r: 



G( P ) = 






P 21 
-8/3 



I. A 4+3( \45nJ P 

For I = 3/4, to = 2, we finally have 

16 / 45tt \3/ 4 3/2 




G(p) = 



25 

9 - ( 
25 \45vr 



(<¥■)'■ v 

3 

V 



9 ( 8 \4/3 n _ 8/3 n> 




(28) 



(29) 



(30) 



We now assume that K is similar to G, but with a 
different value of to. Good results were obtained with to = 
7, so that the linear torque decreases slower towards lower 
viscosities than the horseshoe drag increases. This leads to 
a maximum corotation torque that is in fact larger than 
the unsaturated horseshoe drag, a feature that is already 
apparent in the numerical results shown in Fig. [2] (see also 



3.0 [ 



2.5 - 




0.0 L 
0.1 



Figure 3. Corotation torque as a function of viscosity (through 
the p-parameter). Symbols indicate numerical results, obtained 
with FARGO, for a q = 1.26 • 10 planet in a locally isothermal 
disc with a = 1/2, f) = 1, h = 0.05 and b/h = 0.4. The solid 
curve denotes the model of equation |32| 



Baruteau & Lin ( 2010 1). For this different value of to, we 
find 



K(p) 



1- 



16 (457r\3/4 3/2 
25 I 28 I V 

( 28 \ 4 /3, 



P 



-8/3 




(31) 



We can now construct the complete model for the coro- 
tation torque, including both saturation and the cut-off at 
high viscosity: 



r c ,baro = G(p)F(p)rhs,baro + (1 — -R'(p))Fc,lin,baro , 



(32) 



which uses equations |22[ |30| and |31[ and the torques from 
equations [4] and [6] We compare the complete model to nu- 
merical results in Fig.|3j where we have used the expressions 
for the horseshoe drag and linear corotation torque from Pa- 
per I. As in the isothermal case, for the numerical results, 
x s was determined using a dichotomic search of the sepa- 
ratrices, with a; s being calculated as the geometric average 
of the horseshoe half- widths at ±1 rad (Casoli & Masset 



20091. Equation 32 captures the main features of the nu- 



merical results quite well. We have checked that this is true 
for different density and temperature profiles. 



5 NON-ISOTHERMAL DISCS 

5.1 Disc models 

We now consider the behaviour of non-barotropic effects un- 
der the influence of diffusion. The disc models are similar to 
those in the isothermal case, but in addition we solve the 
energy equation. We stress again that we do not include 
viscous energy dissipation in the equations; we assume that 
dissipation, together with radiative cooling, sets the back- 
ground temperature profile. 
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5.2 Thermal diffusion 

Since we now have an extra physical quantity that is evolved, 
namely the temperature, we consider thermal diffusion in 
addition to viscous diffusion. Since we want the disc to be in 
thermal equilibrium, i.e. viscous heating should compensate 
for radiative losses, we can constrain the range of values we 
need to consider for the thermal diffusion coefficient. 

Radiative diffusion gives rise to a thermal conductivity 



K R = 



acT 6 



4aT' 



(33) 



3np 3np 

where a is the radiation constant, c the speed of light, a 
the Stefan-Boltzmann constant and k the opacity. The cor- 
responding thermal diffusion coefficient \ IS 

Kr _ 4( 7 - l)aT 4 _ 47(7 - f )crT 4 
c v p 

c vP T = P/( 7 - 



X 



3npP 3np 2 H 2 Q 2 ' 

where we have used the internal energy e 

1), with c v the specific heat at constant volume. 



(34) 



Thermal equilibrium requires that (Kley & Crida 2008) 



<rT efl = -fEfi 2 , 



where T e g is the effective temperature (Hubeny|1990[| 



T c s — T /r cf r, 

with an effective optical depth 



3r y/3 

Tcff = y + x 



f 

4V 



with r = KpH is the vertical optical depth. 
Using equation [35] in equation |34| leads to 

X ™ -1 9 , .. /. 2\/3 2 
1/ 4 V 3r 3t^ 



(35) 
(36) 

(37) 
(38) 



where Pr is the Prandtl number, and we have taken E = 
2pH. We therefore have that in the limit of r — > 00, Pr will 
be of order unity, with in general Pr < f . In other words, 
thermal and viscous diffusion will be of equal magnitude 
in the optically thick inner regions of protoplanetary discs. 
In the outer regions, thermal diffusion will come to domi- 
nate, leading eventually to isothermal behaviour for r«l 
( Paardekooper fc Mellema||2006a I . Therefore, we only have 
to consider cases for which \ ^ v i or > equivalently, Pr Sj I. 
Note that this also holds in the presence of other heating 
sources such as external irradiation: in thermal equilibrium, 
any non-viscous heating source will require an increase in \ 
to radiate this energy away. Additional sources of heating 
therefore necessarily decrease Pr even further. 

We first consider a barotropic disc, including thermal 
diffusion, and show that thermal diffusion does not affect 
the barotropic part of the corotation torque. Next, we study 
the case v v = Xpi with \p the thermal diffusion coefficient 
at the location of the planet, and show how this case is 
different to the barotropic case. We then adjust the model 
to incorporate the case v p ^ Xp- 



5.3 Constant entropy disc 

A disc with P = P(E) oc E 7 (or constant entropy) initially 
behaves in a similar way to an isothermal disc. Since the 
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Figure 4. Total torque on a q = 1.26 - 

embedded in a h = 0.05 disc with a 

7 = 7/5, making £ = 0. Two values for thermal diffusion are 

considered: \p = (solid line), and Xp = 10 — 6 fpf2p (dotted 

line). Results were obtained with RODEO. 



entropy is constant, there is no entropy-related horseshoe 
drag, and since vortensity is conserved along streamlines 
in barotropic flows, the vortensity-related horseshoe drag is 
exactly the same as in the isothermal case. The only differ- 
ences between the isothermal and constant entropy cases is 
the value of 7, which affects the total Lindblad and corota- 
tion torques (see Paper I), but also the saturation properties 
through its effect on x s . Since x s decreases for increasing 7, 
for a given value of v p , the horseshoe drag in a constant en- 
tropy disc with 7 > 1 will appear less saturated compared 
to the isothermal case. 

We show in Fig. [1] that the time evolution of the torque 
in a constant entropy disc is not affected by thermal diffu- 
sion. Since we do not apply viscous heating or radiative cool- 
ing, entropy is conserved for \p — 0, keeping the disc, and 
the associated torque on the planet, fully barotropic. When 
applying thermal diffusion, entropy is no longer strictly con- 
served, but since both thermal and viscous diffusion act to 
restore the original density and temperature profiles, the 
disc essentially remains barotropic. The torque on the planet 
is therefore very similar in both cases. This is not true for 
very large values of \p : when thermal diffusion is strong 
enough to affect the horseshoe turn itself, we approach the 
isothermal regime. This limit is considered in section [5~5] 



5.4 The case z/ p = Xp 

We now consider the case Pr = 1. We have used the same 
disc model as for Fig. [3] but now solving the energy equa- 
tion with 7 = 7/5. Note that this disc model has negative 
gradients in both entropy and vortensity, making the to- 
tal horseshoe drag positive, and dominated by the entropy- 
related part (see Paper I). The results are displayed in Fig. 
[5] The solid curve denotes the model of equation|32| with the 
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results were obtained with p x = 3p„-\/xp/^p/2, or 
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Figure 5. Corotation torque as a function of viscosity (through 
the p-parameter). Symbols indicate numerical results for a q = 
1.26 ■ 10~ 5 planet in a disc with a = 1/2, = 1, h = 0.05, 
7 = 7/5 and b/h = 0.4. The Prandtl number equals unity in all 
runs. The solid curve denotes the model of equation |32[ with the 
appropriate values of the barotropic and non-barotropic linear 
and non-linear unsaturated torques, the dashed curve the modi- 
fied model discussed in section[5.4|(see also equation |51|. 



horseshoe drag and linear corotation torque from Paper I. 
It is clear that the torque saturates faster than predicted by 
equation |32| This has to do with the fact that the entropy- 
related part of the horseshoe drag essentially originates from 
a single streamline close to the separatrix. At this location, 
advection is at its maximum efficiency, since the velocity 
with respect to the planet is largest, and diffusion has a hard 
time to keep up. At the same time, the entropy-related part 
needs both thermal and viscous diffusion to remain unsatu- 
rated. Thermal diffusion alone can not provide fresh angular 
momentum to the horseshoe region. Even if thermal diffu- 
sion (or entropy diffusion) can restore the entropy gradient 
sufficiently, if viscosity is not strong enough to diffuse the re- 
sulting vortensity production at the separatrices, the torque 
will still saturate ( Paardekooper fc Papaloizou|2008 l. In this 
sense, the entropy torque can be doubly saturated. 

Rather than try to incorporate these effects into the 
simple isothermal model presented in section |4~T) we try to 
adjust the model, guided by the numerical results to allow 
for the effects described above. Both can be incorporated by 
making the substitution 



F(p) -► F{p v )F{p x 



(39) 



where p x oc Pvy/Xpfvp is the saturation parameter related 
to thermal diffusion and p v is the saturation parameter as- 
sociated with viscosity, in equation |32| Since advection is 
at its most efficient in the region where the entropy-related 
torque is generated, we expect p x > p v for v p — Xp- Good 



Px 



.20 T 3 

piipj. s 

2ttx p 



(40) 



We can see indeed that the dashed curve in Fig. [5] obtained 
with equation |39| gives a good fit to the data at large p v . 
Note that the substitution equation [39] needs to be made for 
the entropy-related part of the horseshoe drag only. 

The cut-off at high v and x should also be modified. 
Note that in this regime, advection being at its most efficient 
actually helps to keep the torque non-linear at higher levels 
of diffusion. Rather than defining new functions G and K, 
we have found that the numerical results can be fitted by 
making the substitutions 



G( P ) -»■ >/G(p,)G(Px) 

and 

1 - K(p) -* ^(1-K(p„))(l-K(p x )). 



(41) 



(42) 



Here again, these substitutions only concern the entropy- 
related part of the horseshoe drag. The resulting model is 
shown by the dashed curve in Fig. [5] Note that at higher lev- 
els of diffusion, the RODEO results differ by 20% from results 
obtained with FARGO. We have verified that half of this is 
due to the different diffusion implementation: entropy diffu- 
sion implemented in RODEO gave lower values for the torque 
by approximately 10% at high diffusion. The remaining 10% 
can be attributed to a higher value of the unsaturated torque 
using RODEO (see Paper I). This dependence on the nature 
of the diffusion should be kept in mind when interpreting 
these results in terms of radiative discs, where x is due to 
radiative diffusion. The error in the model compared to the 
simulations can thus be up to 20 %. 



5.5 The general case z/ p $C Xp 

We now consider the general case, still with the restriction 
Pr ^ 1. We have run a grid of models with different v p an 
X P , and the results are displayed in Fig. [6] All numerical 
results were obtained with RODEO, using thermal diffusion. 

For all values of p v , Pr = 1 gives the lowest value for 
the corotation torque. Going up in thermal diffusion, at a 
fixed value of p v , the torque increases towards the original 
model, denoted by the solid curve. This is consistent with 
the effect discussed in the previous section, that a higher 
thermal diffusion coefficient is needed to fight advection of 
entropy along the separatrix. For all values of Xp> the im- 
proved model gives a reasonable fit to the numerical results, 
with all deviations less than approximately 20 %. 

For Xp > 10 -5 , the corotation torque starts to decrease 
again. This is due to thermal diffusion starting to affect the 
horseshoe turns. For Xp large enough, we expect to enter 
the locally isothermal regime, where the corotation torque 
is given by the barotropic, non-linear horseshoe drag, plus 
the linear, entropy-related corotation torque. Since the lin- 
ear corotation torque is in general smaller than its non-linear 
counterpart, the corotation torque will decrease towards 
higher values of Xp- This is illustrated in Fig. W\ where we 
show the time evolution of the corotation torque at a fixed 
value of Up — 10 for different values of Xp- The torque 
decreases with Xp, but the barotropic part of the corotation 
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Figure 6. Corotation torque as a function of viscosity (through the p-parameter). Symbols indicate numerical results for a q = 1.26 • 10 -5 
planet in a disc with a = 1/2, = 1, h = 0.05, 7 = 7/5 and b/h = 0.4, for different values of thermal diffusion \ p . The solid curve 
denotes the model of equation |32| (saturation of the non-barotropic horseshoe drag is similar to the barotropic part), with the appropriate 
values for the barotropic and non-barotropic parts of the linear and non-linear unsaturated torques. The remaining curves indicate the 
modified model discussed in section [5. 4| for different values of Xp ( see a l so equation [5TI . Different panels highlight different values of Xp> 
with the corresponding model prediction. Numerical results were obtained with RODEO. 



torque remains non linear. The barotropic horseshoe drag 
plus the linear entropy-related corotation torque amount to 
1.4IV However, at these high values of Xp, we also need to 
consider the effect of thermal diffusion on the effective value 
of 7, which will approach unity for large values of Xp- This 
we discuss next. 



5.6 Lindblad torque 

The Lindblad torque is due to a superposition of linear 
waves, and does not strongly depend on viscosity, for rea- 
sonable values of the viscosity v. Thermal diffusion, on the 
other hand, can play an important role, and is necessary to 
recover the isothermal result for large values of the thermal 
diffusivity Xp- The difference between the isothermal and 
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for different values of thermal diffusion Xp- Results were obtained 
with RODEO. 
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adiabatic Lindblad torque basically results from a difference 
in sound speeds. Since Fl oc c~ 2 , rL,i so /rL,adi = 7- 

Consider the dispersion relation for plane waves oc 
exp(i(cut - kx)) ( |Mihalas fc Mihalas|1984) : 

ixk 2 /co 



2,2 1 



(43) 



1 — ijxk 2 / u ' 
with c s the adiabatic sound speed. From this, one can derive 



2Q 



= ,(44) 



7 Q + \^J{i 2 Q 2 - l) 2 + (4Q - 2 7 Q) 2 + 2 7 2 Q 

with Q = x<^/c 2 . In the limit of no thermal diffusion (Q = 
0), v p — > c s : waves propagate at the adiabatic sound speed. 
For Q — > 00, v p — > Cs/"/: waves propagate at the isothermal 
sound speed. 

The torque for any Fourier component will now scale 
as Vp 2 , which makes the isothermal torque a factor of 7 
stronger than the adiabatic torque. In principle, ui = mfi p 
for the problem at hand, and one should therefore consider 
the effects of diffusion on each Fourier component separately. 
For simplicity, we consider a single, representative value of m 
to describe the effects of diffusion. Since most of the torque 
comes from m ~ 2/3h, we set 



2Xp^p 
3hc' 2 



2Xp 



3h a r 2 n v 



It is convenient to define an effective 7: 



7cff 



2Q7 



(45) 



(46) 



lQ + \^2^/( n 2 Q 2 + l) 2 - 16Q 2 (t^I) + 2 7 2 Q 2 ~ 
so that 

7cffr L /r = -2.5- 1.7/3 + 0.1q, (47) 

which then replaces equation[3] Note that the horseshoe drag 
also scales with "/ e g, since the horseshoe width depends on 
7 ( Paardekooper fc Papaloizou|2009b \ . Equation 46 is plot- 
ted in Fig. 18] as a function of a x = Xp/(c b H), for three 
different scale heights h. It is clear that a substantial ther- 
mal diffusion coefficient is required in order to modify the 
Lindblad torque (typically Xp — 10~ 3 r 2 f2 p ). We have over- 
plotted numerical results obtained with RODEO, for a disc 
with h = 0.05, a = 3/2, /3 = 1 and 7 = 5/3. Since this 
disc has no entropy or vortensity gradients, the corotation 
torque vanishes. We can therefore measure 7 e g from the total 
torque, which should be given by equation [3] with 7 = 7 c ft. 
Since for equation [46] we have chosen a single value of m 
to signal the transition from isothermal to adiabatic, we ex- 
pect the numerical results to show a smoother transition 
between 7 ff = 5/3 and 7 c ff = 1. This is indeed what is 
observed in Fig. [8] Note, however, that the model does a 
good job in predicting where the transition takes place. We 
have checked that this also holds for different values of h. 
Around this transition, the maximum error in 7 e a observed 
is 20 %. Note that, since all unsaturated torques are propor- 
tional to 7^ (see below), the error in j e g only affects the 
magnitude of the total torque, not the sign. The dependence 
of 7eff on thermal diffusion will therefore not affect planet 
migration in a qualitative way, but it is necessary to connect 
the adiabatic to the isothermal regime. 



5.7 A modified torque formula 

Before we combine the results of the previous sections 
into a single torque formula, we restate the assumptions 
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made in deriving these results. We have worked in the two- 
dimensional approximation, considered laminar discs only, 
and kept the planet on a fixed circular orbit. We have ig- 
nored the self-gravity of the disc, assumed thermal equilib- 
rium, and considered a simplified energy equation, ft is clear 
that future studies should be aimed at clarifying the role of 
magnetic turbulence and radiative energy transport in keep- 
ing the corotation torque unsaturated. 

We have mainly focused on a g = 1.26- 1CP 5 using a soft- 
ening parameter b/h = 0.4, which seems to be an appropri- 
ate value to mimic 3D effects (Paardekooper & Papaloizou 



2009a|). In the model, the only place where the mass and 



the softening parameter enter is in the values of the unsat- 
urated torques, which have a q 2 dependence and depend on 



b/h (Paardekooper & Papaloizou 2009a I, and in the value 



of x s , which can be written as (Paardekooper & Papaloizou 



2009b) 



= C{b/h) s 



(48) 



where C can be approximated by a power law around b/h = 
0.4: 



c 



1.1 



0.4 
b/h 



-1/4 



(49) 



(see Paper I). The dependence on softening of the linear 
torques was studied in Paardekooper & Papaloizou ( 2009a I , 



and of x B and the non-barotropic horseshoe drag in Paper I. 
Given a different value of the softening, it is straightforward 
to adjust the formulae below by using equations |48| and |49| 
The value of x B then features in the saturation parameters 
p v (given by equation 19 1 and p x (equation 40 1. Note in 



particular that the functional form of F, G and K will not 
change for different values of b/h. 

The total torque on a low-mass planet is given by the 
Lindblad torque plus the corotation torque: 



r L + r c 



(50) 



with Fl given by equation[3] but with 7 — > ■y c g (see equation 
47}. 

The corotation torque is split into a barotropic part and 
an entropy-related part: 



r c = r r 



TO + Pc 



(51) 



The barotropic part of the horseshoe drag is not affected by 
thermal diffusion, therefore 



= FhB,b a roF(p„)G(p,y) + (1 — K(p v ))T ^lin.ba 



(52) 



where rhs,baro is the fully unsaturated horseshoe drag (equa- 
tionB] but with 7 — > j e g), F c ,im,baro is the linear barotropic 
corotation torque (equation [6] b ut with 7 — > 7 c ft), F (p) gov- 
erns saturation (see equationJ22J| , and G(p) and K(p) govern 
the cut-off at high viscosity (see equations 30 and 31 1. 

For the non-barotropic, entropy-related corotation 
torque, we need the modifications discussed in section [5~4| 

r c , c „t = r hs , cnt F(p„)F( Px )VG(p,,)G(p x ) + 



V(l - K(pu))(l - K(p x ))r c ,nn,ent , 



(53) 



where Fh s ,cnt is the fully unsaturated, entropy-related part of 
the horseshoe drag (equation[5j but with 7 — > 7 c ft), r C: ii n>cnt 
is the linear, entropy-related part of the corotation torque 



(equation V7\ but with 7 — > 7 e g), and p x is the saturation 
parameter associated with thermal diffusion (see equation 
401. 

From Fig. [61 we estimate that the maximum error in 
the corotation torque is approximately 20 % compared to 
the simulations. From Paper I, we know that the error in 
the Lindblad torque is smaller in general, so for the total 
torque we expect differences of 20 % between model and 
simulations as well, except close to the transition between 
isothermal and adiabatic behaviour, where the error in the 
Lindblad torque can be 20 % as well (see Fig. [8}. 



6 DISCUSSION 

We have presented numerical results and some simple mod- 
els that capture the behaviour of the torque due to Lindblad 
resonances and horseshoe drag on low-mass planets embed- 
ded in gaseous discs in the presence of viscous and thermal 
diffusion. Several simplifications of the real physical system 
were made to keep the problem tractable. 

First of all, we have worked with vertical integrated 
quantities only, which makes the problem two-dimensional. 
The entropy-related horseshoe drag seems to play an impor- 
tant role in three-dimensional discs as well ( |Kley et al.|2 009), 
and it has been shown that, in an isothermal set-up, the ef- 
fect of the horseshoe drag is stronger in three dimensions 
(Masset et al. 2006). Clearly, a three-dimensional model of 
the horseshoe region is desirable. This would also get rid 
of the dependence of the torque on the softening parameter 
that is necessarily part of a two-dimensional model (see Pa- 
per I). The results of Kley et al. (20091 indicate, however, 



that the main features of three-dimensional simulations can 
be reproduced nicely within the two-dimensional approxi- 
mation. 

We have only considered laminar, viscous discs. In re- 
ality, discs are expected to be turbulent in regions where 



the gas couples to the magnetic field of the disc (Balbus 
fc Hawley]|1991 1. The interaction of turbulence with horse- 



shoe turns is an area that is still largely unexplored. For the 
isothermal case, with a simplified turbulence model based 



on stochastic forcing, Baruteau & Lin (20101 showed that 



turbulent models show similar behaviour as laminar mod- 
els. It remains to be seen how full magnetic turbulence, in a 
non-isothermal setting, affects the horseshoe drag. 

We have modelled radiative cooling as thermal diffu- 
sion, with a constant thermal diffusion coefficient. In ra- 
diative models, \ wm depend strongly on temperature 
and opacity. Moreover, radiative diffusion will probably be 
most efficient in the vertical direction, while in our two- 
dimensional approximation thermal diffusion is restricted to 
the plane of the orbit of the planet. It was noted in |Kley fc| 
Crida ( 2008 1 that this is in fact an unfavourable case for sus- 



taining the corotation torque compared to vertical cooling 
in radiative models. Note that we have chosen \ such that 
the thermal diffusion time scale is identical to the vertical 
cooling time scale. Since both act to restore the original tem- 
perature profile, it is to be expected that including heating 
and cooling (without thermal diffusion) will lead to simi- 
lar effects as thermal diffusion alone. When acting together, 
they may reduce the effective Prandtl number since the flow 
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then has two ways of restoring the original temperature pro- 
file. 

We have taken the background gradients of tempera- 
ture and surface density to be given. In reality, they will be 
determined by a balance of radiative heating and cooling 
and viscous dissipation. It was shown recently in |Lyra et al.| 
(20101 that discs where these profiles are calculated self- 



consistently allow for outward Type I migration to occur for 
a large fraction of the life time of the disc. Of course, since 
the mass of the disc decreases with time, radiative cooling 
becomes more and more efficient, leading in the end to a 
situation where the entropy-related torques are fully linear. 
The strong dependence of the horseshoe drag on the 
local temperature gradient makes it relatively easy to trap 
planets at special locations in the disc. While in the isother- 
mal case, a surface density jump ( Masset et al.|2006 1 or at 



least a positive surface density gradient ( Paardekooper & 



|Papaloizou 2009a) is required, in the non-isothermal case a 
trap is readily established if a local negative temperature 
gradient exists. Even if the disc as a whole does not permit 
outward Type I migration, low-mass planets may stop their 
inward journey when reaching such a temperature trap. 

If the disc parameters allow for outward migration, a 
natural planet trap occurs at the transition between the op- 
tically thick and optically thin regions of the disc. Outside 
this radius, we have isothermal, and therefore inward, Type 
I migration, while inside this radius, we have adiabatic, and 
therefore outward, Type I migration. The evolution of this 
equilibrium radius was explored in more detail in |Lyra et al.| 
(20101. Because of the dependence of the corotation torque 



on the saturation parameter p, which depends on the planet 
mass through x s , different planet masses will have different 
equilibrium radii. The impact of this migration behaviour 
on the gravitational interactions between planets are still 
largely unexplored. 

We have focused on low-mass planets, comparable to 
the Earth. While planets of approximately 1 Jupiter mass 
start to open up deep annular gaps in the disc ( |Lin fe| 
Papaloizou 19861, there exists an interesting intermediate 



regime where planets comparable to Neptune experience an 
enhanced ho rseshoe drag (jMasset et al.|2006J|. Th is was in- 
terpreted in Paardekooper & Papaloizou (2009b I as being 



due to the ineffectiveness of the Lindblad wake to push 
against horseshoe region, leading to a larger value of the 
horseshoe width at these masses. This means that these 
planets are more readily subject to outward migration, lead- 
ing to mass segregation even without considering saturation. 
The planet was kept on a fixed circular orbit through- 
out this paper, and we have ignored any radial viscous flow 
in the gas. If the relative radial speed of the planet with 
respect to the gas is large enough, the streamline topology 
will be modified, leading to a flow structure similar to that 



encountered in Type III migration ( |Masset fe Papaloizou 
2003| |Pepliiiski et al.||2008[ ). Although in this case the dis- 



cussion on the corotation torque will need to be modified, 
we comment that a very massive disc, with a Toomre Q of 
order unity, is needed for a low-mass planet to migrate fast 
enough to achieve this situation. 

Finally, we have ignored the self-gravity of the disc. The 
impact of self-gravity on Type I migration was discussed in 
|Pierens fe Huri| ( [20051 ) and |Baruteau fc Masset| ( |2008b[ ), 
finding an acceleration of inward migration due to a shift in 



the Lindblad resonances. This should be taken into account 
especially in massive discs. A quantitative estimate of this 
effect can be inferred from Fig. 8 of [Baruteau fc Masset] 
[ |2008b| >. 



7 CONCLUSIONS 

We have presented a simple model that can capture the com- 
plicated behaviour of the corotation torque with respect to 
viscous and thermal diffusion. The formulae presented in 
section |5.7| show good agreement with two-dimensional hy- 
drodynamic simulations to within 20 %. The strong depen- 
dence of the total torque on temperature gradients and the 
level of viscous and thermal diffusion allow for a wealth of 
possible outcomes of Type I migration, depending on global 
disc parameters. 
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